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For the calculation of multi-loop Feynman integrals, a novel numerical method, the Direct Com- 
putation Method (DCM) is developed. It is a combination of a numerical integration and a series 
extrapolation. In principle, DCM can handle diagrams of arbitrary internal masses and external 
momenta, and can calculate integrals with general numerator function. As an example of the 
performance of DCM, a numerical computation of two-loop box diagrams is presented. Further 
discussion is given on the choice of control parameters in DCM. This method will be an indis- 
pensable tool for the higher order radiative correction when it is tested for a wider class of physical 
parameters and the separation of divergence is done automatically. 
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1. Introduction 

The high-statistics data in high-energy physics requires the theoretical prediction with enough 
accuracy. The prediction can be given by perturbative calculation in quantum field theory. Then 
the multi-loop integral is an indispensable component for the theoretical study. 

We define the multi-loop integrals by the following formula where the space-time dimension is 
denoted as n = 4-28K Here we confine the discussion to scalar integrals only. Since the method 
presented below is basically numerical, the inclusion of a numerator will be straight-forward. 



r - d n £j JL l 



where the propagator is D r = q 2 .—m 2 + ie, N is the number of propagators and L is the number of 
loops. We combine the propagators by the standard Feynman parameter integral 

and perform integration with respect to the loop momenta. We obtain 

T(N-nL/2) r g(l-Exr) 

{4n) nL/2 X/ ' V J ll dX >- U n/2 {V _ i£) N-nL/2> 

V=M 2 -— , M 2 = Y,x r m 2 r (1.4) 

r 

where U and W are polynomials in the x parameters[|l]]. 

In Section 2, we propose a unique method to calculate the integral / in Eq. 1.3. We call the 
method the Direct Computation Method (DCM). In the preceding works[Q], DCM has successfully 
calculated one-loop and two-loop diagrams. As an example we show the results for the two-loop 
box diagrams in Section 3 and also report a study on the parameters in DCM. We discuss further 
aspects of DCM in Section 4. 

2. Method 

DCM consists of a regularized numerical integration and an extrapolation of a numerical se- 
quence. 

In order to illustrate the idea of regularized numerical integration, let us consider a simple 
integral: 

dx 

m 2 — sx{ \ — x) — is 

As is shown in FigHJ, the integrand of Eq.|2.l| has singular points when s > Am 2 . Analytically, 
this can be handled by taking e as an infinitesimal positive quantity, or, in other words, by the hyper- 
function formula l/(z — ie) = P(l /z) + i7l8 (z). Numerical computation is unstable if the integrand 
is divergent at some points. A way to solve the situation is to deform the path in the complex x 



JO ' 



The symbol e is reserved for the (infinitesimal) parameter in the propagator. 
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Figure 1: Integral paths for (a) analytical calculation and for (b) DCM. The cross (x) stands for the singu- 
larity of integrand. 



plane to avoid the singular points. However, the deformation would be very complicated for the 
integrand of a multi- variable function. Another possibility is to assume that £ is a finite quantity as 
is shown in Fig.[j](b). Then, the numerical integration along the x-axis is stable. If we can calculate 
the limiting value for e — > 0, it is the value of the integral J. 
The loop integral 1(e) is defined by 



/(e) 



if / 



8{\-^x r )(V + ie) N - nL l 2 

Tjn/2(y2 + £ 2}N-nL/2 ' 



(2.2) 



and if £ and 8 are finite, the numerical integration can be performed using a suitable numerical 
computation library. 

We use £ determined by a (scaled) geometric sequence 



£ = £, = £ /(A c )', (/ = 0,1,. 
for constants So,A c (A c > 1). Then we expect 

/ = lim /(£/). 



(2.3) 



(2.4) 



Repeating the numerical integration, we obtain a sequence of numerical values of /(£/) for 
I = 0, 1, . . . Jmax- From these values and using an extrapolation method, we can estimate the value 
of / with enough accuracy. 

Next we discuss the singularity originating from 8 — > 0. In Eq.[2.2[, V 2 + £ 2 is positive and 



U is positive semi-definite. Only at the boundary of the integration region U becomes . When 



8 — > +0, it is either integrable like / dxdy- 

Jo ( 



1 



or non 



■integrable like / dx x . In the latter 

Jo x l 6 



(x+ y y- s 

case it develops a pole term ~ 1 /8 as the ultraviolet singularity 3 . Depending on the masses and 
external momenta, V can be inside the integration region to develop the imaginary part of / 4 , 
and also can be at the boundary of integral region (as in the case of U) to develop an infrared 



2 U is a sum of monomials of x. 

3 The singular pole also appears in the first factor of Eq 



1.3 



if/V-nL/2<Ofor<5 = 0. 



For illustration, one assumes that N — nL/2 = 1 and the variables are transformed into V and x' r variables. Then, 



omitting the Jacobian and other details, the imaginary part becomes 
the inner integral is finite if V2 > > V\ and otherwise. 



dV 



V 2 + e 2 



In the limit as e 



hO, 
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singularity pole ~ 1/5. If Eq ]2.2| is free from these singularities, we just put 5=0. If not, the 
integral in Eq. |2.2| is denoted as 1(8, e) and we first calculate 1(8) = lim/^oo/(5,£/) as in Eq.2.4 
fixing 5. Then, we assume the following form: 



7(5) = .-. + ^-l + C + C 1 5 + --- 
o 



(2.5) 



We calculate 1(8) for a set of values of 5 and estimate the coefficients Cy. For instance, in case of 
a single pole, 5/(5) = C_i + 0)5 + 0(8 2 ) and we extract C_i and Co 5 . 

So much for the description of DCM and one can understand the necessity of an efficient 
library for the numerical integration and that for the extrapolation. For the former we use DQAGE[Q] 
which is a variant of Gaussian quadrature. Since it works adaptively, one can specify the accuracy 
of the numerical results, although the high accuracy costs in computation time. For the latter we 
use Wynn's e algorithm[Q] 6 which predicts the limiting value by the following iteration. We set 
the results of the numerical integration as initial values of the series a(I,k): 



a(/,-l)=0, a(l,0) =/(£,), 1 = 0,1,-- 
The element a(l,k + 1) is obtained by the following recurrence relation: 



a(l,k+\) = a(l + \,k-\) + 



1 



a[l +l,k) — a(l,k) 



1 = 0,1,- 



(2.6) 



(2.7) 



Whilst the a(l,k) values with odd k are meant to store temporary numbers, the a(l,k) with even k 
give extrapolated estimates. 

3. Numerical results 



We calculate the integrals for the two-loop box diagrams shown in Fig.||. The parameters are 
m i = m 2 = m 5 = m 6 = m = 50GeV, «?3 = 1114 = mj = M = 90GeV, p\ = p\ = p\ = p\ = m 2 and 
t = (p\+ p-i) 1 = — (100) 2 GeV 2 . We take s = (p\ + Pi) 2 variable and introduce a dimensionless 
variable f s = s/m 2 . 



(a) 
Pi : 



P2 
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(b) 
Pi : 
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Figure 2: (a)Two-loop planar box diagram and (b) Two-loop non-planar box diagram. The mass of internal 
line k is m^. Each external momentum flows inward. 



The explicit form of the U and W functions is found in [g]. Since there is no ultraviolet/infrared 
divergence, we put 5 = in Eqjl.3|. The integral is 6-dimensional and we perform a transformation 



5 And also C\ is necessary if the first factor of Eq. 1 .3 is singular. 



This '£' has nothing to do with the parameter in propagators. 
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of the integration variables onto a 6-dimensional hypercube [0, l] 6 . By this transformation one can 
cancel common variables between the numerator and the denominator. The results are presented in 
Fig.|| and in Fig.Q In [§J, we verified the results by comparison with another computation which is 
a combination of algebraic transformations and numerical integration, and also by the consistency 
check between the real and imaginary parts through the dispersion relation. 




Figure 3: Numerical results of / for the planar diagram in units of 10~ 12 GeV~ 6 for 0.0 < f s < 25.0 and 
t = — 10000. OGeV 2 . Plotted points are the real part (bullets) and the imaginary part (squares). For the latter, 
the region f s > 10 is not yet computed. 




Figure 4: Numerical results of / for the non- planar diagram in units of 10~ 12 GeV" 6 for 0.0 < f s < 20.0 
and t = — 10000. OGeV 2 . Plotted points are the real part (bullets) and the imaginary part (squares). 

In DCM, the validity of the extrapolation depends on the choice of the values of e. We keep 
finite value of e in the numerical integration, so that its physical dimension is the same as the 
squared mass. We have two parameters in Eq. [2.3[ A c is normally A c = 2 and we can use A c = 1 .2 
or 1.3 to obtain a less computational intensive sequence of integrals as e decreases more slowly. In 
order to study the choice of Eq, the initial value of the iteration, we have calculated the following 
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0.01 1 100 10000 le+06 

e 

Figure 5: The real part(lower points) and imaginary part(upper points) of one-loop scalar vertex integral 
1(e) are shown. The horizontal red line is the exact (analytical) value. 



on-mass-shell one-loop vertex integral as an example: 



Jo<x+ y <i M 2 (\-x-y)+m 2 e (x + y) 2 -sxy-ie 

Here, s = 500 2 GeV 2 and m e = 0.5 x 10 _3 GeV. The value of this integral is computed with a given 
value of M and e = 1 .2 70 ~ m for m = 0, 1, . . . , 120. Then we use the values for m, m + 1 , . . . ,m + 14 
as the input of extrapolation. This means that we take l max = 14 and £o = 1.2 70_m GeV 2 for m = 
0,1,..., 100. 

Table.l Extrapolated values of /. 15 values, /(£/) = Eq/(1.2) 1 , (I = 0, ...,14), are used for the 
extrapolation. Error is not the difference from the analytical value but estimated from the extrapo- 
lation. 







Real part of / 


error 


Imaginary part of / 


error 


3.49E+05 


-1 


75104242540072E-05 


1.65E- 


09 


2.2555636002093 1E-09 


1.80E- 


09 


5.63E+04 


-1 


75105247961553E-05 


3.02E- 


13 


-2.96789310051 170E-05 


8.78E- 


08 


9.10E+03 


-1 


75105248407057E-05 


1.55E- 


14 


4.35402513918612E-05 


1.82E- 


09 


1.47E+03 


-1 


751 052509934 12E-05 


6.78E- 


21 


4.34982757936633E-05 


4.04E- 


13 


2.37E+02 


-1 


75105250996915E-05 


2.76E- 


16 


4.34982757936633E-05 


4.04E- 


13 


3.83E+01 


-1 


75105250989527E-05 


1.40E- 


17 


4.34982788205288E-05 


1.92E- 


17 


6.19E+00 


-1 


75105250987944E-05 


3.58E- 


18 


4.34982788206820E-05 


5.55E- 


16 


1.00E+00 


-1 


75105251 1041 17E-05 


1.22E- 


15 


4.34982788654878E-05 


9.93E- 


14 


1.62E-01 


-1 


75105250921691E-05 


1.98E- 


16 


4.34982793676347E-05 


2.1 IE- 


14 


2.61E-02 


-1 


75105251352866E-05 


1.71E- 


15 


4.34982790401757E-05 


4.97E- 


14 


4.21E-03 


-1 


75 10525 1110498E-05 


4.07E- 


15 


4.34982788582295E-05 


2.21E- 


17 


(analytical) 


-1.75105250974494E-05 






4.34982788194091E-05 







The calculated results for M = 90GeV are shown in Fig.|| and Table.l. It is to be noted that e 
in DCM is obviously finite. One can see that even when the value of 1(e) differs from the analytical 
value, the extrapolation gives a good estimation. There is a finite region that shows the agreement 
between the analytical value and the extrapolated one. This behavior demonstrates that DCM is 
stable up to some extent for the choice of Sq parameters. We have tested the similar analysis for 
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several values of M (M = 1, 10 1 , 10 2 , 10 3 GeV) and found similar behavior. Though we need more 
tests for this point, we can temporally conclude that £o is not need to be very small but it can be set 
to the typical squared mass in the denominator of the integrand. 

In the calculation of the two-loop box diagrams described above, we have checked the conver- 
gence of the extrapolation step-by-step. The value used for £o is 1.2 40 ~ 1.2 45 GeV 2 which would 
be consistent with the above conjecture. 

4. Summary 

In this paper, we have outlined DCM and calculated the two-loop scalar box integral as an 
example to show its applicability. Since the radiative correction in the electroweak theory (or in the 
SUSY model) involves various combinations of mass parameters in the integrand, DCM is a good 
candidate to handle general loop integrals. 

In order to use DCM for the calculation of higher-order radiative corrections, we plan to pro- 
ceed to the following research. 

1. The method should be tested for a wider class of diagrams with various combinations of 
masses and external momenta and with the numerator structure. 

2. Further study on the choice of parameters Eo,A c is required for a stable application. 

3. The variable transformation is sometimes important for good convergence. This is to be 
processed in an automatic manner. 

4. In dimensional regularization, the ultraviolet/infrared divergence appears as a pole 1/5. The 
separation of the infrared pole is already done successfully in A similar treatment of 
ultraviolet poles is expected. 

5. Sometimes DCM needs long computational time. It will be important to perform the com- 
putations in a parallel computing environment. 
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